Assignment #2

library(tidyverse)         # for graphing and data cleaning
library(tidymodels)        # for modeling
library(stacks)            # for stacking models
library(naniar)            # for analyzing missing values
library(lubridate)         # for date manipulation
library(moderndive)        # for King County housing data
library(vip)               # for variable importance plots
library(DALEX)             # for model interpretation
library(DALEXtra)          # for extension of DALEX
library(patchwork)         # for combining plots nicely

data("lending_club")

Put it on GitHub!

Here is my GitHub link.

Modeling

We’ll be using the lending_club dataset from the modeldata library, which is part of tidymodels. The outcome we are interested in predicting is Class. And according to the dataset’s help page, its values are “either ‘good’ (meaning that the loan was fully paid back or currently on-time) or ‘bad’ (charged off, defaulted, of 21-120 days late)”.

Tasks

1. Explore the data, concentrating on examining distributions of variables and examining missing values.

lending_club %>% 
  ggplot(aes(x = funded_amnt)) +
  geom_density() + 
  facet_wrap(vars(Class))

lending_club %>% 
  ggplot(aes(x = int_rate)) +
  geom_density() + 
  facet_wrap(vars(Class))

lending_club %>% 
  ggplot(aes(x = annual_inc)) +
  geom_density() + 
  facet_wrap(vars(Class))

lending_club %>% 
  ggplot(aes(x = addr_state, fill = Class)) +
  geom_bar(position = "fill")

lending_club %>% 
  count(Class)

lending_club %>% 
  group_by(addr_state) %>% 
  summarize(count = n()) %>% 
  arrange(desc(count))

2. Do any data cleaning steps that need to happen before the model is build. For example, you might remove any variables that mean the same thing as the response variable (not sure if that happens here), get rid of rows where all variables have missing values, etc.

create_more_bad <- lending_club %>% 
  filter(Class == "bad") %>% 
  sample_n(size = 3000, replace = TRUE)

lending_club_mod <- lending_club %>% 
  bind_rows(create_more_bad)

3. Split the data into training and test, putting 75% in the training data.

set.seed(494)

lending_split <- initial_split(lending_club_mod,
                             prop = 0.75)
lending_train <- training(lending_split)
lending_test  <- testing(lending_split)

4. Set up the recipe and the pre-processing steps to build a lasso model. Some steps you should take:

• Make all integer variables numeric (I’d highly recommend using step_mutate_at() or this will be a lot of code). We’ll want to do this for the model interpretation we’ll do later.
• Think about grouping factor variables with many levels.
• Make categorical variables dummy variables (make sure NOT to do this to the outcome variable).
• Normalize quantitative variables.

lending_recipe <- recipe(Class ~ .,
                         data = lending_train) %>% 
  step_rm(acc_now_delinq, delinq_amnt) %>% 
  step_mutate_at(all_numeric(),
                 fn = ~as.numeric(.)) %>% 
  # step_mutate(annual_inc =
  #             case_when(annual_inc <= 9875                          ~ 10,
  #                       annual_inc > 9875   && annual_inc <= 40125  ~ 12,
  #                       annual_inc > 40125  && annual_inc <= 85525  ~ 22,
  #                       annual_inc > 85525  && annual_inc <= 163300 ~ 24,
  #                       annual_inc > 163300 && annual_inc <= 207350 ~ 32,
  #                       annual_inc > 207350 && annual_inc <= 518400 ~ 35,
  #                       annual_inc > 518400                         ~ 37)) %>%
  # step_mutate(annual_inc = as.factor(annual_inc)) %>% 
  step_normalize(all_predictors(), -all_nominal()) %>% 
  step_dummy(all_nominal(), -all_outcomes())

lending_recipe %>% 
  prep(lending_train) %>% 
  juice()

5. Set up the LASSO model and workflow. We will tune the penalty parameter.

lasso_lending_mod <- logistic_reg(mixture = 1) %>%
  set_engine("glmnet") %>% 
    set_args(penalty = tune()) %>% 
      set_mode("classification")

lasso_lending_mod

## Logistic Regression Model Specification (classification)
## 
## Main Arguments:
##   penalty = tune()
##   mixture = 1
## 
## Computational engine: glmnet

lending_workflow <- workflow() %>%
  add_recipe(lending_recipe) %>%
    add_model(lasso_lending_mod)

lending_workflow

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: logistic_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 4 Recipe Steps
## 
## ● step_rm()
## ● step_mutate_at()
## ● step_normalize()
## ● step_dummy()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Logistic Regression Model Specification (classification)
## 
## Main Arguments:
##   penalty = tune()
##   mixture = 1
## 
## Computational engine: glmnet

6. Set up the model tuning for the penalty parameter. Be sure to add the control_stack_grid() for the control argument so we can use these results later when we stack. Find the accuracy and area under the roc curve for the model with the best tuning parameter. Use 5-fold cv.

set.seed(494)

lending_cv <- vfold_cv(lending_train, v = 5)

lending_lasso_pen_grid <- grid_regular(penalty(), levels = 10)

ctrl_grid <- control_stack_grid()

lending_lasso_tune <- lending_workflow %>% 
  tune_grid(resamples = lending_cv,
            grid = lending_lasso_pen_grid,
            control = ctrl_grid)

lending_lasso_tune %>% 
  show_best(metric = "accuracy")

lending_lasso_tune %>% 
  collect_metrics() %>% 
  filter(.config == "Preprocessor1_Model01")

7. Set up the recipe and the pre-processing steps to build a random forest model. You shouldn’t have to do as many steps. The only step you should need to do is making all integers numeric.

ranger_recipe <- recipe(Class ~ .,
                         data = lending_train) %>% 
  # step_rm(acc_now_delinq, delinq_amnt) %>% 
  step_mutate_at(all_numeric(),
                 fn = ~as.numeric(.))
  # step_mutate(annual_inc =
  #             case_when(annual_inc <= 9875                          ~ 10,
  #                       annual_inc > 9875   && annual_inc <= 40125  ~ 12,
  #                       annual_inc > 40125  && annual_inc <= 85525  ~ 22,
  #                       annual_inc > 85525  && annual_inc <= 163300 ~ 24,
  #                       annual_inc > 163300 && annual_inc <= 207350 ~ 32,
  #                       annual_inc > 207350 && annual_inc <= 518400 ~ 35,
  #                       annual_inc > 518400                         ~ 37)) %>%
  #step_mutate(annual_inc = as.factor(annual_inc)

8. Set up the random forest model and workflow. We will tune the mtry and min_n parameters and set the number of trees, trees, to 100 (otherwise the next steps take too long).

ranger_spec <- rand_forest(mtry = tune(),
                           min_n = tune(),
                           trees = 100) %>% 
  set_mode("classification") %>% 
  set_engine("ranger")

ranger_spec

## Random Forest Model Specification (classification)
## 
## Main Arguments:
##   mtry = tune()
##   trees = 100
##   min_n = tune()
## 
## Computational engine: ranger

ranger_workflow <- workflow() %>% 
  add_recipe(ranger_recipe) %>% 
  add_model(ranger_spec)

ranger_workflow

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: rand_forest()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 1 Recipe Step
## 
## ● step_mutate_at()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Random Forest Model Specification (classification)
## 
## Main Arguments:
##   mtry = tune()
##   trees = 100
##   min_n = tune()
## 
## Computational engine: ranger

9. Set up the model tuning for both the mtry and min_n parameters. Be sure to add the control_stack_grid() for the control argument so we can use these results later when we stack. Use only 3 levels in the grid. For the mtry parameter, you need to put finalize(mtry(), lending_training %>% select(-Class)) in as an argument instead of just mtry(), where lending_training is the name of your training data. This is because the mtry() grid will otherwise have unknowns in it. This part can take a while to run.

set.seed(494)

lending_rf_grid <- grid_regular(min_n(), 
                                finalize(mtry(), lending_train %>% select(-Class)),
                                levels = 3) 

ctrl_res <- control_stack_grid()

ranger_cv <- ranger_workflow %>% 
  tune_grid(resamples = lending_cv,
            grid = lending_rf_grid,
            control = ctrl_res)

10. Find the best tuning parameters. What is the are the accuracy and area under the ROC curve for the model with those tuning parameters?

ranger_cv %>% 
  show_best(metric = "accuracy")

ranger_cv %>% 
  collect_metrics() %>% 
  filter(.config == "Preprocessor1_Model4")

The best tuning parameters are 11 for mtry and 2 for min_n.

11. Use functions from the DALEX and DALEXtra libraries to create a histogram and boxplot of the residuals from the training data. How do they look? Any interesting behavior?

best_param <- lending_lasso_tune %>% 
  select_best(metric = "accuracy")
lending_lasso_final_wf <- lending_workflow %>% 
  finalize_workflow(best_param)
lending_lasso_final_mod <- lending_lasso_final_wf %>% 
  fit(data = lending_train)

best_param_ranger <- ranger_cv %>% 
  select_best(metric = "accuracy")

lending_ranger_final_wf <- ranger_workflow %>% 
  finalize_workflow(best_param_ranger)

set.seed(494)
ranger_fit <- lending_ranger_final_wf %>% 
  fit(lending_train)

lending_lasso_explain <- explain_tidymodels(
  model = lending_lasso_final_mod,
  data = lending_train %>% select(-Class),
  y = as.numeric(lending_train$Class == "good"),
  label = "LASSO",
  type = "classification"
)

## Preparation of a new explainer is initiated
##   -> model label       :  LASSO 
##   -> data              :  9643  rows  22  cols 
##   -> data              :  tibble converted into a data.frame 
##   -> target variable   :  9643  values 
##   -> predict function  :  yhat.workflow  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package tidymodels , ver. 0.1.2 , task classification (  default  ) 
##   -> model_info        :  type set to  classification 
##   -> predicted values  :  numerical, min =  0.02883008 , mean =  0.7276827 , max =  0.9883398  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -0.9244458 , mean =  -4.632243e-06 , max =  0.9711699  
##   A new explainer has been created! 

lending_rf_explain <- explain_tidymodels(
  model = ranger_fit,
  data = lending_train %>% select(-Class),
  y = as.numeric(lending_train$Class == "good"),
  label = "Random Forest",
  type = "classification"
)

## Preparation of a new explainer is initiated
##   -> model label       :  Random Forest 
##   -> data              :  9643  rows  22  cols 
##   -> data              :  tibble converted into a data.frame 
##   -> target variable   :  9643  values 
##   -> predict function  :  yhat.workflow  will be used (  default  )
##   -> predicted values  :  No value for predict function target column. (  default  )
##   -> model_info        :  package tidymodels , ver. 0.1.2 , task classification (  default  ) 
##   -> model_info        :  type set to  classification 
##   -> predicted values  :  numerical, min =  0 , mean =  0.7038608 , max =  1  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -0.4 , mean =  0.02381728 , max =  0.29  
##   A new explainer has been created! 

lending_lasso_model_perf <- model_performance(lending_lasso_explain)
lending_lasso_model_perf

## Measures for:  classification
## recall     : 0.9381502 
## precision  : 0.7804386 
## f1         : 0.852058 
## accuracy   : 0.7629368 
## auc        : 0.7772458
## 
## Residuals:
##          0%         10%         20%         30%         40%         50% 
## -0.92444585 -0.68117176 -0.48463606  0.04695045  0.09244145  0.12651880 
##         60%         70%         80%         90%        100% 
##  0.16297193  0.20979491  0.29607744  0.41138573  0.97116992

lending_rf_model_perf <- model_performance(lending_rf_explain)
lending_rf_model_perf

## Measures for:  classification
## recall     : 1 
## precision  : 1 
## f1         : 1 
## accuracy   : 1 
## auc        : 1
## 
## Residuals:
##    0%   10%   20%   30%   40%   50%   60%   70%   80%   90%  100% 
## -0.40 -0.01  0.00  0.00  0.00  0.01  0.02  0.03  0.05  0.08  0.29

hist_plot <- 
  plot(#lending_lasso_model_perf,
       lending_rf_model_perf, 
       geom = "histogram")
box_plot <-
  plot(#lending_lasso_model_perf,
       lending_rf_model_perf, 
       geom = "boxplot")

hist_plot + box_plot

The histogram is fairly normally distributed with a slight right skew towards positive residuals. The boxplot shows that there is a large amount of residual variance which is evidence by the RMSE being located outside of the IQR, as well as the long tail on the right.

12. Use DALEX functions to create a variable importance plot from this model. What are the most important variables?

set.seed(494)

lending_lasso_var_imp <- 
  model_parts(
    lending_lasso_explain
    )
lend_lasso_VIP_plot <- 
  plot(lending_lasso_var_imp, show_boxplots = TRUE)

lending_rf_var_imp <- 
  model_parts(
    lending_rf_explain
    )
lend_rf_VIP_plot <- 
  plot(lending_rf_var_imp, show_boxplots = TRUE)

lend_lasso_VIP_plot

lend_rf_VIP_plot

Interest rate seems to be the most important variable by a large margin for the LASSO. For the random forest, annual income and interest rate are also much more important than any of the other variables, with annual income being more important than interest rate, which is the opposite of the LASSO.

13. Write a function called cp_profile to make a CP profile. The function will take an explainer, a new observation, and a variable name as its arguments and create a CP profile for a quantitative predictor variable. You will need to use the predict_profile() function inside the function you create - put the variable name there so the plotting part is easier. You’ll also want to use aes_string() rather than aes() and quote the variables. Use the cp_profile() function to create one CP profile of your choosing. Be sure to choose a variable that is numeric, not integer. There seem to be issues with those that I’m looking into.

cp_profile <- function(explainer, newObs, varName) {
  predict_profile(explainer = explainer,
                  new_observation = newObs,
                  variables = varName) %>% 
    rename(yhat = `_yhat_`) %>% 
    ggplot(aes_string(x = varName, y = 'yhat')) + 
    geom_point()
}

set.seed(494)

cp_profile(explainer = lending_rf_explain, 
           newObs = lending_train %>% slice_sample(), 
           varName = 'annual_inc')

14. Use DALEX functions to create partial dependence plots (with the CP profiles in gray) for the 3-4 most important variables. If the important variables are categorical, you can instead make a CP profile for 3 observations in the dataset and discuss how you could go about constructing a partial dependence plot for a categorical variable (you don’t have to code it, but you can if you want an extra challenge). If it ever gives you an error that says, “Error: Can’t convert from VARIABLE to VARIABLE due to loss of precision”, then remove that variable from the list. I seem to have figured out why it’s doing that, but I don’t know how to fix it yet.

rf_pdp <- model_profile(explainer = lending_rf_explain, 
                        variables = c('int_rate', 'annual_inc', 'open_il_12m', 'open_il_24m'))

plot(rf_pdp,
     variables = 'int_rate',
     geom = 'profiles')

plot(rf_pdp,
     variables = 'annual_inc',
     geom = 'profiles')

plot(rf_pdp,
     variables = 'open_il_12m',
     geom = 'profiles')

plot(rf_pdp,
     variables = 'open_il_24m',
     geom = 'profiles')

15. Fit one more model type of your choosing that will feed into the stacking model.

lending_knn_mod <-
  nearest_neighbor(
    neighbors = tune("k")
  ) %>%
  set_engine("kknn") %>% 
  set_mode("classification")

lending_knn_workflow <- workflow() %>% 
  add_model(lending_knn_mod) %>% 
  add_recipe(lending_recipe)

lending_knn_tune <- lending_knn_workflow %>% 
  tune_grid(resamples = lending_cv, 
            grid = 4,
            control = ctrl_grid)

16. Create a model stack with the candidate models from the previous parts of the exercise and use the blend_predictions() function to find the coefficients of the stacked model. Create a plot examining the performance metrics for the different penalty parameters to assure you have captured the best one. If not, adjust the penalty. (HINT: use the autoplot() function). Which models are contributing most?

set.seed(494)

lending_stack <-
  stacks() %>% 
  add_candidates(ranger_cv) %>% 
  add_candidates(lending_lasso_tune) %>% 
  add_candidates(lending_knn_tune)

## Warning: Predictions from the candidates c(".pred_bad_lending_lasso_tune_1_02",
## ".pred_bad_lending_lasso_tune_1_03", ".pred_bad_lending_lasso_tune_1_04",
## ".pred_bad_lending_lasso_tune_1_05", ".pred_bad_lending_lasso_tune_1_06",
## ".pred_good_lending_lasso_tune_1_02", ".pred_good_lending_lasso_tune_1_03",
## ".pred_good_lending_lasso_tune_1_04", ".pred_good_lending_lasso_tune_1_05",
## ".pred_good_lending_lasso_tune_1_06") were identical to those from existing
## candidates and were removed from the data stack.

lending_blend <- lending_stack %>% 
  blend_predictions()

lending_blend

## ── A stacked ensemble model ─────────────────────────────────────

## 
## Out of 18 possible candidate members, the ensemble retained 3.
## Lasso penalty: 0.001.

## 
## The 3 highest weighted member classes are:

## # A tibble: 3 x 3
##   member                          type             weight
##   <chr>                           <chr>             <dbl>
## 1 .pred_good_ranger_cv_1_4        rand_forest      12.7  
## 2 .pred_good_lending_knn_tune_1_1 nearest_neighbor  0.930
## 3 .pred_good_ranger_cv_1_7        rand_forest       0.297

## 
## Members have not yet been fitted with `fit_members()`.

autoplot(lending_blend)

autoplot(lending_blend, type = "weights")

The random forest models are contributing the most, with the knn model contributing slightly.

17. Fit the final stacked model using fit_members(). Apply the model to the test data and report the accuracy and area under the curve. Create a graph of the ROC and construct a confusion matrix. Comment on what you see. Save this final model using the saveRDS() function - see the Use the model section of the tidymodels intro. We are going to use the model in the next part. You’ll want to save it in the folder where you create your shiny app.

lending_final_stack <- lending_blend %>% 
  fit_members()

lending_final_stack

## ── A stacked ensemble model ─────────────────────────────────────

## 
## Out of 18 possible candidate members, the ensemble retained 3.
## Lasso penalty: 0.001.

## 
## The 3 highest weighted member classes are:

## # A tibble: 3 x 3
##   member                          type             weight
##   <chr>                           <chr>             <dbl>
## 1 .pred_good_ranger_cv_1_4        rand_forest      12.7  
## 2 .pred_good_lending_knn_tune_1_1 nearest_neighbor  0.930
## 3 .pred_good_ranger_cv_1_7        rand_forest       0.297

saveRDS(lending_final_stack, file = "lending_final_stack")

lending_stack_test <- lending_test %>% 
  bind_cols(predict(lending_final_stack, new_data = lending_test, type = "prob")) %>% 
  bind_cols(predict(lending_final_stack, new_data = lending_test))

lending_stack_test %>% 
  accuracy(.pred_class, Class)

lending_stack_test %>% 
  roc_auc(Class, .pred_bad)

autoplot(roc_curve(lending_stack_test, Class, .pred_bad))

Shiny App

The shiny app is not finished but a link to the GitHub page is provided to show progress.

• Here is the GitHub repository for my Shiny app.

Coded Bias

• The part of the film that impacted me the most was the part where the high school aged kid was stopped on the street in London because the facial recognition technology misidentified him due to his race. It provided me evidence for how much of a problem this truly is. I was most surprised by how widespread facial recognition was used in China, and how China uses it to create a citizen rating. I was also surprised at how the United States does this without nearly as much transparency, which is a huge problem in my eyes. It was truly shocking how transparent China is about facial recognition software uses, and how some citizens are okay with what China is doing. I was also shocked that there is no AI or Machine Learning Algorithm regulation in the United States. Overall, I felt generally unsettled by the documentary, as it made me fear for the future of algorithm usage in the United States.